Impurity effects at finite temperature 
in the two-dimensional 5 = 1/2 Heisenberg antiferromagnet 



Kaj H. Hoglund 1 and Anders W. Sandvik 1,2 

'Department of Physics, Abo Akademi University, Porthansgatan 3, FIN-20500, Turku, Finland 
^Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, U.S.A. 

(Dated: February 2, 2008) 

We discuss effects of various impurities on the magnetic susceptibility and the specific heat of the 
quantum 5=1/2 Heisenberg antiferromagnet on a two-dimensional square lattice. For impurities 
with spin St > (here Si = 1/2 in the case of a vacancy or an added spin, and Si = 1 for a 
spin coupled ferromagnetically to its neighbors), our quantum Monte Carlo simulations confirm 
a classical-like Curie susceptibility contribution S? /3T, which originates from an alignment of the 
impurity spin with the local Neel order. In addition, we find a logarithmically divergent contribution, 
which we attribute to fluctuations transverse to the local Neel vector. We also study frustrated and 
nonfrustrated bond impurities with Si = 0. For a simple intuitive picture of the impurity problem, 
we discuss an effective few-spin model that can distinguish between the different impurities and 
reproduces the leading-order simulation data over a wide temperature range. 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Cx, 75.40.Mg 



I. INTRODUCTION 

The problem of impurities in low-dimensional quan- 
tum antiferromagnets has attracted considerable atten- 
tion ever since the discovery of high-temperature super- 
conductivity in the cupratesi At low concentration, holes 
doped into the Cu02 planes are localized, or have very 
low mobility, and hence static impurities are relevant 
for understanding the initial reduction of antiferromag- 
netism upon doping compounds such as La2Cu04. These 
impurities are expected to be magnetically frustrated^ 
Although not directly related to the breakdown of antifer- 
romagnetism associated with the onset of superconduc- 
tivity, static nonfrustrating impurities, e.g., inert holes 
corresponding to substitution of Cu atoms by nonmag- 
netic Znj&4 also can give important information pertain- 
ing to the nature of the interactions in the CuC>2 planes. 
The same applies to related cuprates where the planes are 
broken up into chains^ or ladders^ Very recently, similar 
impurity problems were also suggested to be of relevance 
to possible physical realizations of quantum computers^ 

On the theoretical side, Heisenberg impurity mod- 
els can be studied by a wide range of modern quan- 
tum many-body methods. Importantly, numerical tech- 
niques, such as quantum Monte Carlo and the density 
matrix renormalization group, can give approximation- 
free results against which analytical approaches can be 
tested on a quantitative level. Once such a program 
has been completed, the applicability of a Heisenberg de- 
scription to an experimental system can be judged with- 
out concerns about approximations in the calculations. 
There is already a large body of work devoted to var- 
ious impurity effects, and a coherent picture is emerg- 
ing. Restricting ourselves to work on single impurities, 
we note several ground state calculations for Heisenberg 
chains ^LULIL ladders jiSiiiLii and the two-dimensional 
(2D) square latticei 13 i 15 i 16 i 17 i 18 i 19 i 20 i 21 i 22 Extensive work 



on chaino 23 : 24 ! 25 and 2D systemo 21 : 26 : 27 : 28 ! 29 at finite tem- 
perature has also been carried out. In this paper, we 
continue our studies 27 of finite-temperature effects of iso- 
lated static impurities in the standard 2D Heisenberg 
model, and also present some results for the 3D system. 

In a recent comprehensive quantum field-theoretical 
workjSi a low-temperature theory of an arbitrary quan- 
tum impurity in a 2D antiferromagnetic host system was 
developed, with the host being either in the T = 
magnetically ordered phase (i.e., in the renormalized- 
classical regime at T > 0) or close to a quantum-critical 
point. In the magnetically ordered phase a leading-order 
classical-like Curie contribution to the impurity suscep- 
tibility was predicted to stem from the coupling of the 
impurity moment Si to the local Neel order as the tem- 
perature T —>■ 0, i.e., xf mp — > Sf/3T. Stimulated in part 
by these theoretical predictions by Sachdev et a2.*2L we 
recently carried out a large-scale quantum Monte Carlo 
(QMC) studjiS 7 . of the 2D S = 1/2 Heisenberg anti- 
ferromagnet and confirmed the Curie prefactor 1/12 in 
the renormalized-classical regime, for a vacancy (miss- 
ing spin) as well as for an added Si = 1/2 impurity 
spin. However, we also discovered a low-T logarithmi- 
cally divergent subleading contribution to the impurity 
susceptibility. This anomaly was attributed to the trans- 
verse component, for which a T-independent behavior 
had been predictedi 21 : 30 Related logarithmic divergences 
had also previously been found, e.g., in an exact study 
of an impurity in the classical 2D Heisenberg model^ 
and in a Green's function treatment of the 2D quantum 
model with an extra spin at T = In the latter study, 
the frequency dependent transverse impurity susceptibil- 
ity XimpC^ 1 = 0, <*;) was found to be log divergent when 
lu — ► 0. More recently, an anomalous susceptibility was 
also found in the Heisenberg model with a finite impurity 
concentration^ 

Recent efforts, by Sachdev and Vojta 2 ^ and Sushkov} 2 ^ 
to explain our previous numerical findings^ have re- 
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suited in more complete field-theoretical descriptions 
where the impurity susceptibility indeed acquires a pre- 
viously unnoticed subleading log divergent contribution. 
Its principal cause can also in these analytical treatments 
be considered as associated with the transverse compo- 
nent, although there are also higher-order logarithmic 
contributions arising from longitudinal fluctuations^2*2il 
In the formulation by Sachdev and Vojta^ which relies 
on an expansion of the nonlinear a model around dimen- 
sionality d = 1, a very detailed form for the logarithmic 
corrections to the impurity susceptibility was given. For 
general impurity spin S^. 
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where p s is the spin stiffness of the bulk-ordered anti- 
ferromagnet in the absence of impurities and the un- 
known constants are in general nonuniversal, but 
become universal when a quantum critical point is ap- 
proached. The first subleading term oc ln(l/T) in Eq. Q 
hence accounts for the log divergent behavior observed in 
our numerical studies. The quantitative agreement be- 
tween Eq. Q and our numerical data is quite remark- 
able, as will be shown in this paper. Our results also 
agree qualitatively with the analytical results obtained 
by SushkoviSS 

The purpose of this paper is to give a more complete 
numerical account of the effects of different types of sin- 
gle static impurities on the magnetic susceptibility of the 
2D 5=1/2 Heisenberg antiferromagnet on a square lat- 
tice. Some of the results were previously summarized 
in Ref. I27I The impurity effects were there determined 
for a vacancy and an added-spin impurity, by calculating 
impurity susceptibilities with the stochastic series expan- 
sion (SSE) QMC techniquei 33 i 34 The impurity suscepti- 
bility is simply the difference between the susceptibilities 
of the pure and doped Heisenberg models. In this paper 
we compare our numerical results for the vacancy and 
added-spin impurity models with the theoretical expres- 
sion in Eq. Q). We also consider an impurity consisting 
of a spin coupled ferromagnetically to its four nearest 
neighbors [see Fig. ^d)]. This coupling arrangement is 
nonfrustrating and can be expected to lead to an Si = 1 
impurity, in contrast to the Si — 1/2 vacancy and added- 
spin impurities. It was suggested by Aharony et aim 
that hole doping the parent compounds of the cuprate 
superconductors could lead to effective frustrated ferro- 
magnetic exchange couplings between nearest neighbor 
Cu spins. Motivated by this scenario, we have also con- 
sidered an impurity model with a single ferromagnetic 
bond, and compared this with a missing bond. The 2D 
Heisenberg antiferromagnet with two vacancies on dif- 
ferent sublattices, and at different separations, is also 
studied in order to further elucidate the behavior of the 
single- vacancy impurity susceptibility. Finally, we have 



considered a single vacancy in the three-dimensional (3D) 
Heisenberg antiferromagnet, for which analytical limiting 
expressions has also been obtained recently. 28 Although 
the main focus of this paper is on the susceptibility, we 
will also present some results for impurity effects on the 
internal energy and the specific heat. 

In Ref. |53 we also introduced effective models for the 
vacancy and added-spin systems. These models are very 
simple few-spin systems constructed in order to capture 
the leading-order impurity effects — they do not contain 
the log corrections. They provide simple physical pic- 
tures of the dominant mechanisms at play in the full 
Heisenberg impurity models. In this paper the effective 
models are discussed in detail, and the concept is further 
demonstrated by results for added-spin impurities with 
different couplings to the host and the ferromagnetically 
coupled in-plane impurity spin. 

The rest of the paper is organized as follows. In Scc.lTTl 
the full Heisenberg and effective models, as well as their 
impurity susceptibilities, are defined. In Sec. lIIII the SSE 
Monte Carlo method is briefly outlined, and the improved 
estimators needed to achieve sufficient statistical accu- 
racy are discussed. We also describe an averaging trick 
used to alleviate the sign problem in our study of the 
frustrated ferromagnetic-bond impurity. The SSE results 
are presented and compared with the corresponding ef- 
fective models in Sec. IIVI Concluding remarks are given 
in Sec.0 In the Appendix we discuss the specific heat of 
the pure Heisenberg model, for which we have obtained 
low-temperature results of unprecedented accuracy. In 
order to provide benchmark results for other calculations, 
we also list some selected high-precision numerical values 
for energies and susceptibilities of systems with different 
impurities. 



II. IMPURITY MODELS AND 
SUSCEPTIBILITIES 

Following Ref. l2ll an impurity susceptibility is defined 
as the difference between the susceptibility of an impurity 
system and the pure system, i.e., 
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where i = b,c,d,e, and /, correspond to the different 
impurity systems shown in Figs.^andEl an d x\ a ) is the 
susceptibility of the pure system. The susceptibilities on 
the right-hand side of Eq. J2J are not normalized by the 
system size, i.e., 
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where the sum is over all the spins of the pure or impurity 
systems. The impurity susceptibilities are hence inten- 
sive differences of extensive quantities, and they provide a 
natural framework for quantifying the effects of different 
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isolated impurities on the susceptibility of the pure sys- 
tem. They also give the leading (linear) dependence on 
the concentration of impurities. The definition in Eq. iffll 
will be used both in the context of the full Heisenberg 
models and the corresponding effective models, both of 
which will be defined in this section. Quantities analo- 
gous to Eq. ||2J) will also be used for the internal energy 
and the specific heat. 

The impurity susceptibility can be separated in com- 
ponents parallel and perpendicular to a given direction. 
Here the separation is done with respect to an axis along 
the orientation of the local Neel order at the impurity. In 
the isotropic 2D Heisenberg antiferromagnet, true long- 
range Neel order sets in, i.e., the spin-rotation symmetry 
of an infinite system is broken, only at T = The 

components xlmp* an d Xhap > where || and _L refer to di- 
rections parallel and perpendicular to the Neel order, are, 
therefore, true physical observables only at T = 0. How- 
ever, our calculations show a temperature behavior that 
confirms an approximate, but conceptually useful, sep- 
aration of the impurity susceptibility in components al- 
ready at low finite T, as will be shown in Sec. IIVI In the 
3D Heisenberg antiferromagnet, Neel order is present al- 
ready at finite temperature, below T c /J ss 0.95j2A which 
makes the two components truly distinguishable. The 
effective impurity models are defined to include a fluctu- 
ating direction given by a classical vector N, describing 
a local Neel order with respect to which susceptibility 
components can be defined. Comparisons with the QMC 
results show that the separation into components is use- 
ful even at a quantitative level. 

A. Full Heisenberg models 

The basis for this study is the isotropic S = 1/2 Heisen- 
berg antiferromagnet on a periodic L x L lattice. This 
model is defined by the Hamiltonian 

N b 

H (") = J ^2 S *(b) ' S J(6)' ( 4 ) 
6=1 

where J > 0, bond b connects the nearest- neighbor sites 
[i(b), j(b)], and Nb is the total number of bonds. Hu) is 
given a pictorial representation to the left in Fig. u3[a), 
and will hereafter be referred to as the full Hamiltonian 
of the pure system. Impurity models are obtained by 
introducing single defects in the pure model. 

We begin by presenting the models with impurity mo- 
ments Si 0. They are illustrated in Fig. 2] When a 
single spin So is removed from the square lattice, the va- 
cancy model shown to the left in Fig. ^b) is obtained. 
We will study it in 2D as well as in 3D. The added-spin 
model, shown to the left in Fig.^c), is obtained by cou- 
pling a single off-plane spin-^ S a antiferromagnetically 
to a spin So on the square lattice. Two different values 
on the coupling strength J± = J and J± = J/2 will be 
considered here. In the limit Jj_ — > oo, the magnetic 
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FIG. 1: Full Heisenberg models and corresponding effective 
models of the (a) pure, (b) vacancy, (c) added-spin, and (d) 
four ferromagnetic bonds systems. Thick solid lines symbolize 
ferromagnetic spin-spin couplings — Jf < 0, with Jf = J- 
The free parameters of the effective models are the couplings 
a and r. 




(f) 



FIG. 2: Full Heisenberg models of the (e) frustrating ferro- 
magnetic bond (Jf = J) and (f) removed bond (Jf = 0) 
systems. 



properties of the added-spin model become equivalent 
to the vacancy model, since the two spins S a and So 
are then locked in a singlet state. An Si = 1 impurity 
is obtained by considering a configuration of four ferro- 
magnetic bonds with one spin in common, as shown in 
Fig. Hid). 

We also consider the models with impurity moments 
Si = shown in Fig. [21 A system with one frustrat- 
ing ferromagnetic bond, where the spin-spin coupling is 
— Jf < 0, is shown in Fig.^e). In the limit Jf/ J — > oo, 
one might hence expect a corresponding Si = 1 impurity 
moment, as the two spins connected by the ferromag- 
netic bond then form a triplet. On the other hand, for 
Jf = clearly Si = 0, and hence a transition between 
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Si = and 1 might be expected at some intermediate 
J p. However, it has been shown that in a Neel ordered 
bulk, the correlations between the spins connected by the 
ferromagnetic bond remain antiferromagnetic^ 7 - which is 
possible because of the broken degeneracy of the triplet 
when it is coupled to an asymmetric environment. Hence, 
counterintuitively, an Si = 1 behavior when T — > may 
actually never be realized with a ferromagnetic bond im- 
purity. We will here consider only the coupling strength 
Jp — J, for which the sign problem due to frustration 
can be alleviated by a position averaging procedure, as 
discussed in Sec. III. Our QMC data shows that the im- 
purity moment Si — in this case, and the behavior is 
similar to the removed-bond impurity model (Jf = 0) 
shown in Fig.E^f). 

Two-impurity models are useful for further clarifying 
the properties of the single-impurity models. Here we will 
consider only the case of two vacancies, which are chosen 
to be either at a fixed short distance from each other or 
maximally separated on the L x L square lattices. 



B. Effective models 

The purpose of introducing effective models is to cap- 
ture the dominant mechanisms at play in the full Heisen- 
berg models with very simple systems containing a mini- 
mal number of adjustable parameters. The effective mod- 
els are constrained by two criteria: (i) they should repro- 
duce the high-T impurity susceptibility, the sign of which 
depends on whether a spin has been added or removed, 
and (ii) they should mimic the expected^ leading-order 
behavior Sf/3T of the impurity susceptibilities at low T, 
i.e., the alignment of the impurity moment with the local 
Neel order. Effective models are here considered for the 
full Si impurity models shown in Fig. ^ 

In our effective models the local Neel order at the im- 
purity is modeled by a classical "nonmagnetic" vector N. 
When a single spin So is removed from the full model of 
the pure system, as shown to the left in Fig.^b), the re- 
maining system has an S — 1/2 ground state due to the 
sublattice asymmetry. Hence, the effective model corre- 
sponding to the vacancy system, shown to the right in 
Fig. |T^ h). is simply defined with a single effective rem- 
nant "environment" spin-i S e . This spin is coupled to a 
classical unit vector N representing the orientation of the 
local Neel order. The magnitude of this order is absorbed 
in the coupling strength r. The effective model for the 
pure system is naturally obtained by reinserting the spin 
So, as shown to the right in Fig. da), with a > for 
antiferromagnetic coupling. The effective model for the 
added-spin system, shown to the right in Fig. ^c) , is ob- 
tained by coupling an extra spin-^ S a to So. Finally, the 
effective model for the system with a configuration of four 
ferromagnetic bonds, shown to the right in Fig. ^d), is 
obtained by simply changing the sign of a, i.e., by making 
the coupling ferromagnetic instead of antiferromagnetic. 
To summarize, the Hamiltonians of the effective models, 



corresponding to the full models in Fig. ^ are 
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rN-S e + aSo • S e , (5a) 

rN • S e) (5b) 

rN • S e + aS ■ S e + J ± S a • S , (5c) 

rN-S e -aS -S e . (5d) 



The parameters r > and a > cannot be derived in 
any trivial way. The magnitude of r should in principle 
depend on T, but the T dependence can be expected to 
be weak once the amplitude of the order has developed 
locally close to the impurity. One could also argue that a 
direct coupling between N and the central spin So should 
be included. Such a coupling is clearly mediated through 
the four nearest neighbors of So- However, in the spirit of 
keeping the models as simple as possible, we here chose to 
accomplish this coupling indirectly through the remnant 
environment spin S e . One can further anticipate that the 
optimum values for the couplings r and a will depend 
on the impurity type, since the effective impurity spin is 
spread out and its coupling is mediated through the local 
environment of So, which will be distorted in different 
ways by different impurities. However, we will show that 
the same parameters, r/J w 1.90 and a/ J ~ 2. 25, 37 
actually give an overall reasonable agreement for all the 
Si > impurity types considered here. 

The procedure for determining the susceptibilities of 
the effective models is straightforward. An external ap- 
plied field h = h z e z defines the z direction. The mag- 
netization operators M(j), corresponding to the effec- 
tive Hamiltonians in Eqs. ©, have the z components 

S? + S* and 



M[ a) = S$ + 
formula 



S!, Mfo = SI 



Mfa = SI 



The susceptibilities are given by the usual 
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where i — a,b,c, and d. Here the inner brackets (•) in- 
dicate the quantum mechanical expectation value for a 
fixed direction of N, and (-)n denotes the classical ori- 
entation average. The imaginary-time evolved operator 



M ( z i} (r) = exp(r J ff ( e ^)M ( 2 i) exp(-riJ ( e i f |). Expressing M { 



(i)J 



in the coordinate system defined by N, 

Mf j} = oostejMjj -sin (e)Mfo , 
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where || and _L denote the directions parallel and per- 
pendicular to N, the expectation values can be easily 
calculated. The second term in Eq. © vanishes. The 
first term can be separated to components: 
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FIG. 3: The components of the susceptibilities of the effective 
models for the pure (a) and vacancy (b) systems. The com- 
ponents of the impurity susceptibility Ximp ~ X(b) ~ X( a ) are 
shown in (c). The inset shows the 4Txjmt? ~ 1/3 behavior as 
T^0. 



where the prefactors originate from the classical orienta- 
tion averaging. The susceptibility components are 

4) = / 1/r ^< M (i)WAfjl ) (0)> = i((Mjl ) ) a ) ) (9a) 

X® = I ' dr(M^(r)M^(0)), (9b) 
Jo 

which can be easily evaluated in the || basis. In this basis 
Eq. i|9a|) has the simple form because [Hf^, A4"JL] = 0. 

A dominant feature of the effective models is the align- 
ment of a quantum spin with a classical vector. This can 
be appreciated by examining the simple effective Hamil- 
tonian HfR, given in Eq. J5bj, for the vacancy system. 



The corresponding susceptibility is given by 

1 || 2 ± 

X(b) - 3%(6) + 7}X(b) 

11 2 1 / r \ 

= Slr + a^fe)- (10) 

The temperature dependence of the two components is 
graphed in Fig. Efb) for r = 1.90. Since S e = l/2 in this 

model, the || component can be written as (1/3)%IL = 

S e 2 /3T. Hence, the classical Curie prefactor S e 2 , instead 
of the usual quantum mechanical prefactor S e (S e + 1), is 
a consequence of the finite coupling between the spin S e 
and the classical vector N. This is precisely the low- 
T leading order behavior proposed^! in Eq. in the 
renormalized classical regime of a 2D antiferromagnet. 
The perpendicular component in Eq. I|ll)|l tends to a con- 
stant at low T. On the other hand, in the limit r — > 0, 
i.e., when S e and N decouple, S e recovers its quantum 
identity and the susceptibility has the usual Curie form 
X z (b) ^S e (S e + l)/3T. 

The effective models also serve the purpose of elucidat- 
ing the steps in determining the impurity susceptibility 
Ximp = Xn,) — Xfay The separation in components of 
the susceptibility for the effective pure system, x\ a y is 
shown in Fig. Efa). The parallel component (l/3)xL) 
vanishes at low T since the spins So and S e are then 
aligned antiferromagnetically with respect to N. The 
perpendicular component (2/3)%^ assumes a constant 
value at low T. The components of the susceptibility 
Xf^-s for the vacancy system were given analytically in 
Eq. (|1L)|) and are shown in Fig.^b). Finally, the compo- 
nents of the impurity susceptibility are shown in Fig.^c). 
At high T, the impurity susceptibility is just the sum of 
the Curie contributions of each independent spin, i.e., 
xfmp } -> 1/4T - 2 /AT = -1/4T. At low T the parallel 
component diverges, while the perpendicular component 
becomes a constant. The inset verifies that the parallel 
component is responsible for the S' e 2 /3T behavior, since 

4^xfrop ~ 1/3 as T — > 0. Besides being capable of repro- 
ducing the expected low-T leading order behavior of the 
full models, the effective models also account quite accu- 
rately for impurity specific behavior at intermediate T, 
as will be shown in Sec. II VI There we also demonstrate 
that the effective models account accurately for the T 
dependence of the internal energy. 



III. QUANTUM MONTE CARLO METHOD 

The numerical method employed here for the full 
Heisenberg models is the operator-loop formulation of 
the stochastic series expansion (SSE) QMC method. It 
is a method based on importance sampling of the terms 
of the Taylor series for the density matrix. Its applica- 
tion to the Heisenberg model has been discussed in detail, 
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e.g., in Refs.|33| and l34l The method is only briefly out- 
lined here in order to discuss some important aspects of 
the impurity work, including a trick for alleviating the 
sign problem for the frustrated-bond system and the use 
of improved estimators for reducing statistical errors. 

The Hamiltonian of the Heisenberg antiferromagnet in 
Eq. (0} can be cast into the form 

T N » TAT 
ff (o) =--^ ir %)|^, (11) 

6=1 

where the operators H i & and i?2,6, defined by 



H- 



'2.6 



: S i(b) S j(b) +S t (b) S t(b)> 



(12a) 
(12b) 



are diagonal and off-diagonal, respectively, in the basis 
{\a) — \Sf, S%, ■ ■ ■ , Sjv)} used in the simulation. An ex- 
act expression for the partition function Z is obtained by 
expanding the density matrix e~@ H in a Taylor series at 
inverse temperature j3 — 1/T (ks — 1). The series can 
be truncated at some expansion power n max = M, since 
terms of order greater than n tx Nf3 give an exponen- 
tially vanishing contribution,™ The truncated partition 
function is then given by 



M 



a,i,bi 



(13) 



Since the matrix element of the operator product takes 
the values and 1, the statistical weight of a contributing 
configuration isS 



W{a,S M ) = 



(-1)" 2 (/3J)"(M -n)\ 
2™M! 



(14) 



A number of M—n identity operators i?o,o = I have been 
inserted in the matrix element of each term in Eq. (|13fl . 
with expansion order n < M, and the change in prefactor 
reflects the number of different ways to distribute the 
n Hamiltonian operators among the M positions. The 
symbol Sm denotes a sequence of operator indices, 



Sm = (ai,h), (a 2 , 62), ■ ■ • , («m, &m), 



(15) 



where a.; 6 {1,2} and b t £ {1, ...,iV&}, corresponding 
to the operators H ai ^ bi in Eqs. lfT2)l . or (a^^) = (0,0), 
corresponding to the identity operator i?o,o- For a given 
sequence Sm the order n then denotes the number of 
non-(0, 0) operators in the sequence. For a nonfrustrated 
lattice, the number 712 of off-diagonal operators (2, bi) in 
the sequence Sm is always even for nonvanishing contri- 
butions, thus yielding a positive definite statistical weight 
W(a,S M ) in Eq. flj. 

With a positive definite expansion, the partition func- 
tion Z can be stochastically evaluated by importance- 
sampling in the configuration space (a, Sm)- For this 



purpose an algorithm consisting of two different config- 
uration updates is used. In the first update (diagonal 
update) the sequence Sm is traversed from beginning 
to end, while attempting substitutions (0,0) <-> (1, bi). 
The substitution (0,0) — * (1, 6») is attempted only if the 
spins connected by bond bi are antiparallel [for a nonva- 
nishing contribution with the definition of the diagonal 
operator in Eq. (|12afl ]. The probabilities to use for ac- 
cepting/rejecting the change have been given elsewhere, 
e.g., in Ref.l34L An accepted attempt changes the expan- 
sion order n by ±1. If an off-diagonal operator (2, bi) is 
encountered no single operator substitution can be car- 
ried out, and instead the saved state \a) is updated by 
flipping the two spins connected by the bond bi, so that 
the state on which the diagonal operators act are always 
available when attempting and update (0, 0) — > (1, 6»). In 
the second update (operator-loop update) the sequence 
Sm is uniquely decomposed into a number Ni of oper- 
ator loops, in which substitutions (1, bi) <-> (2, bi) can 
be carried out, independently with probability 1/2 for 
each loop. All the spins associated with the loops are 
also flipped. During the operator-loop update the order 
n is kept fixed and the weight of the configuration is un- 
changed. The operator- loop update was introduced and 
discussed in detail in Ref. 134 . 

The simulation is started with a random state |a) and 
an empty sequence Sm = (0, 0), (0, 0), ... , (0, 0) of arbi- 
trary (short) length M. One Monte Carlo step (MC step) 
consists of a diagonal update followed by an operator- 
loop update. During the equilibration stage of the simu- 
lation the cutoff M is adjusted to always exceed the max- 
imum order n reached. Hence, the truncated partition 
function Z in Eq. (|13H is no approximation. Observables 
are measured after every MC step and expectation values 
and their errors are determined by the usual method of 
data binning. Estimators for various observables of the 
Heisenberg antiferromagnet, in the context of the SSE 
method, are discussed in Ref. |3jJ The susceptibility is 
given in Eq. where the sum is evaluated in the stored 
state \a). The internal energy and the specific heat are 
given by2£ 



E = - 



C 



(n) 



(16a) 
(16b) 



The operator-loop formulation of the SSE method, as 
described above, is directly applicable to the isotropic 
Heisenberg antiferromagnet. Impurities in the form of va- 
cancies, added spins, and missing bonds can be included 
with only very minor changes in the algorithm. In the 
added-spin impurity case, the only change in a program 
for the pure model is that the acceptance probabilities 
in a diagonal update (0, 0) — > (1, bk), involving the addi- 
tional bond k connecting the impurity spin, depend on 
the bond strength However, the impurities consist- 
ing of a frustrating ferromagnetic bond or four nonfrus- 
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trating ferromagnetic bonds necessitate some additional 
considerations, as will be discussed next. 

For a ferromagnetic bond, the diagonal bond operator 
(|T2a|) is denned as 2(1/4 + SfS*) and the off-diagonal 
(|12b|l is multiplied by —1. During the diagonal update 
the substitution (0,0) — * (l,2>i), where bi is a ferromag- 
netic (antiferromagnetic) bond, is hence attempted only 
if the spins connected by bond bi are parallel (antiparal- 
lcl) . The rules for constructing the operator loops are also 
modified, as discussed in Ref. |3Ja. For the impurity con- 
sisting of four nonfrustrating ferromagnetic bonds, the 
expression for the statistical weight W in Eq. (|14|> is still 
valid if 712 is replaced by the number nA2 of off-diagonal 
operators (2, bi) acting on antiferromagnetic bonds. Be- 
cause of the symmetry of the arrangement of four ferro- 
magnetic bonds, this number also has to be even, and, 
therefore, the weight W is still positive definite. 

A single frustrating ferromagnetic bond (here with cou- 
pling Jf = J) in the Heisenberg antiferromagnet gives 
rise to a sign problem. Proceeding as in the case of four 
ferromagnetic bonds discussed above, the sign would be 
determined by the number of spin flips of antiferromag- 
netic bonds, which now can be even or odd. Since the 
total number of flips still has to be even, we can also 
define the sign as (— 1)™ F2 , where npi is the number of 
spin flips on the ferromagnetic bond. However, we can 
also proceed in a different way which allows for an alle- 
viation of the sign problem by position-averaging when 
Jf — J ■ We then treat the ferromagnetic bond in the 
same way as an antiferromagnetic bond in the diagonal 
update, i.e., a diagonal operator can appear only on an- 
tiparallel spins. The sign will then be given by (— 1)™ F , 
where n f is the total number of operators — diagonal and 
off-diagonal — operating on the ferromagnetic bond. The 
simulation of the system with a ferromagnetic bond then 
proceeds exactly as the simulation of the pure antiferro- 
magnet, i.e., expectation values can be calculated using 
\W\ and by reweighting the measurements with the sign 
S = ( — 1)™ F of the corresponding configuration, 



(A) 



(AS) 



\w\ 



(S) 



\w\ 



(17) 



In practice, however, the calculations become impossible 
when (S)\w\ approaches zero. Here a technique based on 
positional averaging is used to tackle this problem. The 
idea is to replace the sign S of a given configuration with 
the averaged sigri^ 



(18) 



where an average of the sign S(R) = (— l) nF ( R ) is taken 
with respect to all possible locations R of the ferromag- 
netic bond. Expectation values are then given by 
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FIG. 4: SSE results for the magnetic susceptibility of the pure 
2D Heisenberg antiferromagnet for different system sizes L. 
Error bars are smaller than the symbols. The inset shows 
a comparison between the low-T behavior for the pure, va- 
cancy, and added-spin models with L — 64. For the added- 
spin model, two values on the coupling constant J±/J are 
considered. 



This technique was discussed in a more general context 
in Ref. |39j, where it was shown that it significantly al- 
leviates the sign problem of the antiferromagnet with 
randomly positioned ferromagnetic bonds. This came at 
the price of an approximation corresponding to switching 
to an "annealed" disorder. Here, in the single-impurity 
problem, there is no approximation as the trick simply 
corresponds to simultaneously studying systems with all 
possible locations of the ferromagnetic bond. When con- 
sidering only a single position R of the ferromagnetic 
bond, the sign problem will be more severe than with 
a redefinition of the diagonal operator discussed above 
for the four-bond impurity. However, when using the 
position averaging there will be some system size above 
which the statistic is improved. We here obtained expec- 
tation values with reasonable statistical errors for system 
sizes up to L = 32 at temperatures down to T = J/8. 
Since the evaluation of the sign during the simulation is 
completely separate from the sampling procedures, the 
effect of a ferromagnetic bond can actually be obtained 
as a "bonus" while simulating the pure antiferromagnet. 
A drawback of the position averaging method is that it 
does not allow for ferromagnetic bond strengths Jf =/= J, 
except perhaps for Jp very close to J where reweighting 
should work. 

We next briefly comment on the accuracy needed to 
study the impurity effects and the use of improved esti- 
mators for increasing the accuracy. For large L, the effect 
of a single impurity on the magnetic susceptibility \ z is 
very small, as shown in the inset of Fig.0] In order to get 
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acceptable errors for the impurity susceptibilities xf mp in 
Eq. @ very precise values for the individual susceptibil- 
ities are clearly necessary. To achieve this, an improved 
estimator— is used. The general idea is to reduce the 
statistical errors by replacing the value of an observable 
A corresponding to a given Monte Carlo configuration by 
an estimator Ai obtained by averaging over many equal- 
weight configurations during the operator-loop update. 
In the case of the susceptibility, this is particularly sim- 
ple since the magnetization is a conserved quantity. Some 
of the loops will go through (once or multiple times) the 
state | a), i.e., the state on which the ordered operator 
product is acting on in Eq. (|T3|l . Defining of as the sum 
over all the spins in |a) covered by the z:th loop, we 
clearly have 

N Nt 

M* = = £ of. (20) 

We can now average this over all the 2 Nl ways of flipping 
the loops, giving 

x*=^X>*) 2 V ( 21 ) 

Figure 0] shows size-normalized results for the magnetic 
susceptibility obtained this way for the pure 2D Hciscn- 
berg antifcrromagnet, as well as low-T data for systems 
with an impurity. We believe that these results are the 
most accurate ones currently available for this model and 
therefore also list selected numerical data in the Ap- 
pendix. For the internal energy and the specific heat, 
Eqs. I|16a(l and i|16b(l . no improved estimator of the type 
discussed above can be constructed. The energy can nev- 
ertheless be calculated to high accuracy, as also shown in 
the Appendix. For the specific heat, it is very difficult to 
reach good accuracy at low temperature. Nevertheless, 
we are able to clearly discern the expected 41 behavior 
C oc T 2 at low T, as shown in Fig. ^] in the Appendix. 

IV. RESULTS 

Here, in Sec. IIV AI we begin by presenting susceptibil- 
ity results for the Si ^ impurities illustrated in Fig. ^ 
In Sec. IIV Bl we consider the case of a vacancy in a 3D 
system, and in Sec. IIV 01 we look at the system with two 
vacancies. We discuss results for the Si = bond impu- 
rities (Fig. [3J in Sec. HVDl In Sec. IIV El we summarize 
our results for the impurity effects on the energy and the 
specific heat. 

A. Si 7^ impurities 

The impurity susceptibilities for a vacancy and an 
added spin with J± = J are shown in Figs.JSJa) and[SJb), 
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FIG. 5: The impurity susceptibilities of the (a) vacancy and 
(b) added-spin (Jj_ = J) models, for different system sizes 
L. The dotted lines show the expected asymptotic behavior 
4Txf mp — ► 1/3 as T — > 0. Error bars are smaller than the 
symbols. 

respectively. The results are multiplied by 4T. At high 
T the data for different system sizes L coincide, while at 
lower T finite-size effects are clearly seen for L < 16. The 
finite-size effects are due to the S — 1/2 ground states 
of the vacancy and added-spin models, to which the sys- 
tem converges below an L dependent crossover tempera- 
ture, as has recently been discussed by Sushkov^ For the 
largest system size considered here, L — 64, all finite-size 
effects are eliminated within statistical errors for temper- 
atures down to T/J = 1/32. The observed behavior at 
high T for both impurity types is due to the fact that 
the total susceptibility is then just the sum of the Curie 
contributions of each independent spin, i.e., 

z.(b.c) _ z _ z 
Ximp X(6,c) X( a ) 



4T AT AT v ' 

as T — ► oc. The minus (plus) sign is for the vacancy 
(added-spin) impurity model. According to the expres- 
sion in Eq. Q , the leading order behavior of the impu- 
rity susceptibility is ATxf mp ~ ASf/3 as T — > 0. For a 
Si = 1/2 impurity, the constant value 1/3 should then 
be approached at low T. This is also clearly observed 
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FIG. 6: L = 64 results for the impurity susceptibilities of the 
vacancy and added-spin models is compared to the results of 
the corresponding effective models. The dotted line shows the 
value 1/3. 

in the size-converged (L — 64) data for the vacancy im- 
purity, shown in Fig. Ufa) . For the added-spin impurity 

shown in Fig. EJb) , an approach of ^T\l^ to 1/3 is also 
likely, although the convergence occurs at lower T than 
for the vacancy. At intermediate T the results for the two 
different impurity types are strikingly different. Specifi- 
cally, the shoulder-like structure with a minimum around 
T/J w 0.8 observed in the added-spin data has no coun- 
terpart in the vacancy data, but in both cases there is a 
maximum at T/J ss 0.2. Some of the differences clearly 
are related to the different T — > oo behaviors. 

In Fig. the size-converged SSE data are compared 
with the results of the effective models. Results are 
also shown for the added-spin impurity with J± = J/2. 
The values of the two parameters of the effective models, 
a/ J = 2.25 and r/J — 1.90^ were chosen for optimal 
overall agreement between the SSE data and the effective 
model results, for both the vacancy and the added-spin 
systems. For this choice of values, the effective models 
reproduce the added-spin data with a remarkable preci- 
sion down to T/J rs 0.1, for both J± = J and J± = J/2. 
Moreover, with the same set of values a reasonable agree- 
ment is also obtained for the vacancy system. Hence, the 
same parameters describe well a wide range of coupling 
strengths to the added spin (the vacancy corresponds to 
J±/J = oo). 

In each of the three cases shown in Fig. [fjl the effective 
models also reproduce the low-T leading-order behavior 
suggested in Eq. JTJ, i.e., 4Txf mp ~ 4S?/3 = 1/3 for a 
Si = 1/2 impurity. Hence, the effective models clearly 
contain the dominant impurity physics and are able to 
distinguish between different impurity types in a broad T 
range. In analogy with the results for the effective mod- 
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FIG. 7: SSE data for xf mp - 1/12T of the vacancy and added- 
spin models with system sizes L — 64. Straight lines and 
dashed curves are fits of the theoretical results in Ref. [2^ to 
our low-T simulation data. The dotted curve shows the 1/6T 
behavior. 



els, the observed low-T leading-order behavior of the full 
Heisenberg models is ascribed to a susceptibility compo- 
nent parallel to a locally Neel ordered domain coupled to 
the impurity, i.e., (l/3)xL ( p ~ Sf/3T = 1/12T, where 
i = b,c, and Si = 1/2. 

We next examine the thermodynamic low-T impurity 
susceptibilities more closely by subtracting from them 
the leading-order term Sf/3T, The resulting quantities 
should then describe the transverse impurity susceptibil- 
ities at low T, i.e., (2/3)xtp ~ xfn? " S?/3TM The 

results in Fig.0for Ximp — 1/12T of the vacancy impu- 
rity, and of the added-spin impurity with Ji = J, show 
an apparent logarithmically divergent behavior as T — > 0. 
The results for the added-spin impurity with J± = J/2 
are not conclusive in this regard, but a similar log diver- 
gent behavior at still lower temperatures is clearly plau- 
sible. As J±/J — > oo, the magnetic properties of the 
added-spin model should become equivalent to those of 
the vacancy model. In the limit Jj_/J —> 0, on the other 
hand, the added spin is decoupled from its host and the 
impurity susceptibility becomes simply the susceptibility 
of a single spin (1/4T), i.e, xf mp - 1/12T - 1/6T. When 
comparing the SSE results in Fig. with each other, 
it then seems that the log divergent behavior starts at 
higher T as the magnitude of the coupling to the added 
spin, Jj_/J, is increased. This can be naturally under- 
stood as an impurity moment strongly coupled to the 
environment can develop only at T below J±. 

According to the theoretical expression by Sachdev and 
Vojtaj2£ Eq. Q , the slopes of the the low-T curves should 
be equal on the log-linear scale used in Fig. [7J The slope 
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FIG. 8: Impurity susceptibilities for different system sizes L of 
the impurity model with four ferromagnetic bonds. The solid 
curve shows the result of the corresponding effective model, 
which assumes the asymptotic value 4/3 (shown by the dotted 
line) at low T. The log divergent behavior of xfm P — 1/3T, for 
a system of size L = 32, is shown in the inset, where the solid 
line and the dashed curve are fits of the theoretical results in 
Ref. [2^ to our low-T simulation data. 



is given by Sf /3irp s , where Si is the "bare" impurity spin 
and p a is the spin stiffness of the bulk-ordered antifcr- 
romagnet, for which we use the value p 8 jJ = 0.181. 42 
Our results for the vacancy and the Jj_ = J added 
spin are indeed consistent with this prediction. The 
straight solid lines are fits of the leading logarithmic part, 
oc \n(Cip s /T), of Eq. to the low-T data, whereas 
the dashed curves show fits including also the subleading 
correction oc Tln(C2p s /T). For the vacancy system we 
find C\ sa 1.7 (in the leading-order fit) or C\ « 1.6 and 
C 2 ~ 0.3, for the added-spin system [Jy = J) C\ ~ 50 
or C\ ~ 73 and C2 ~ 184. For the added-spin impurity 
with J± = J/2, no fit can be made with only the leading 
term, and we find C\ w 10 5 and C2 ~ 10 19 . 

Results for the impurity model with a configuration of 
four ferromagnetic bonds are shown in Fig. |S1 Again, 
the results for different L coincide at high T, while 
finite-size effects are seen at lower T. The high-T ob- 
served behavior, 4Tx^p — > as T — > 00, is due to 
the fact that the susceptibilities of the doped and the 
pure models cancel, since there is an equal number of in- 
dependent spins in both models at high temperatures. 
The ground state spin of this model is S = 1, and 
hence also an impurity moment Si — 1 can be antici- 
pated. The low-T finite-size susceptibility should then 
be AT X Z ~ 4T[S(S + 1)/3T] = 8/3 and in the thermody- 
namic limit 4Tx z ~ 4T(S' 2 /3T) = 4/3. This behavior is 
indeed seen in Fig. EI f° r L = 4 and 8 the low-T behavior 
dictated by the ground state spin can be observed, while 
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FIG. 9: The impurity susceptibilities of the 3D Heisenberg 
antiferromagnet with a vacancy, for different system sizes L. 
The dotted line shows the value 1/3. A comparison between 
Ximp — 1/12T of the 2D and 3D models is shown in the inset. 



for L = 32 the low-T susceptibility is size-converged at 
least to T/J = 1/16 and is consistent with a convergence 
to 4/3. We also show results for the corresponding ef- 
fective model. Using the same values for a/ J and r/J 
as previously for the vacancy and added-spin effective 
models, changing only the sign of a, the behavior agrees 
qualitatively with the SSE results. The inset of Fig. [S] 
shows xfrnp"* ~ 1/3T, which at low T should be domi- 
nated by the transverse component (2/3)x^p • Again, 
an apparent log divergent trend is observed. The straight 
line and the dashed curve are fits of Eq. QJ. Data for the 
two lowest T are not included in these fit because of the 
finite-size effects that most likely remain here. Never- 
theless, the results support the universal low-T prefactor 
(slope) of the leading logarithmic correction. 



B. Vacancy in a 3D system 

Here we discuss the case of a vacancy in the 3D Heisen- 
berg antiferromagnet. Some predictions — were recently 
made also for this system, but since we have not achieved 
sufficient accuracy they are not tested in detail here. The 
leading-order behavior can nevertheless be extracted. In 
Fig. the SSE data are shown for different system sizes 
L (N = L 3 ), and a comparison between the 3D and 2D 
data is shown in the inset. The high-T behavior, as well 
as the low-T finite-size effects, have the same explana- 
tions as those given for the 2D results. For the largest 
system size, L = 16, most finite-size effects are eliminated 
within statistical errors in the T range considered. The 
observed thermodynamic behavior is reminiscent of the 
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FIG. 10: Impurity susceptibility for different system sizes L 
of the square lattice with two vacancies. The vacancies are 
as far apart as possible in (a). In (b) the L = 32 results 
are shown for the cases when the two vacancies are nearest 
neighbors (open symbols) and at distance r = (2, 1) from each 
other (solid symbols). 



2D results in Fig. |S[a) , with the exception that the tran- 
sition to a constant-valued behavior now occurs abruptly 
at T/J sa 0.95, which is the Neel temperature Tjy of the 
models There are signs of a singular behavior of the im- 
purity susceptibility at the transition. At temperatures 
T < Tn, the susceptibility is seen to follow very closely 
the proposed2ii2& behavior 5 2 /3T. It should be noted 
that although 3D order sets in below T/v, our finite-size 
systems nevertheless do not break the symmetry and the 
direction of the Neel vector is not fixed. In an infinite 
symmetry-broken system, the S 2 /3T behavior would not 
be present if the magnetization fluctuations are defined 
with respect to the average in the direction of the fixed 
Neel vector. 

In the inset of Fig. [51 the perpendicular component 

( 2 / 3 )x,4 P = Xf m p - (l/3)xL p is compared to the analo- 
gous quantity of the 2D model. Although the statistical 
accuracy is not very high at low temperature, it is clear 
that the behaviors are different. The 3D results do not 
indicate any log divergent behavior of the type observed 
in the 2D system. Instead an almost constant behavior 
is observed, as also predicted in the field theory^ 



C. Two vacancies in 2D 

Next we present SSE results for the 2D Heisenberg an- 
tiferromagnet with two vacancies on different sublattices. 
The results for the impurity susceptibility are multiplied 
by a factor 1/2, so that single-impurity values should 
be obtained when the correlation length is much shorter 
than the separation between the vacancies. When T is 
lowered, interactions between the impurities become im- 
portant as the correlation length £ grows exponentially. 
At T corresponding to a correlation length of the same 
order as the vacancy separation, the moments due to the 
two vacancies on different sublattices are pinned by the 
local Neel order antiparallel to each other, resulting in 
a rapid quenching of the parallel component of the im- 
purity susceptibility. Hence, xfmp does not diverge as 
T — ► 0. The data shown in Fig. II Uf a) is for the case of 
maximum separation of two vacancies on different sublat- 
tices; r = (L/2 — l,L/2). Since £ diverges exponentially 
as T — > 0, the point at which xf mp deviates from the di- 
vergent single-vacancy behavior moves only very slowly 
to lower T as L is increased. For larger L, an almost con- 
stant xfmp i s observed. However, no sign of convergence 
of the plateau value is seen. Clearly, in a system of finite 
size there will always be some interaction also between 
the perpendicular components of the two vacancies, and 
hence even for large L the two-vacancy model does not 
trivially reproduce the single- vacancy results below some 
temperature. It is plausible, however, that the roughly 
ln(£) divergence of the plateau height seen in Fig. HUf a) 
continues as L — > oo. This would be fully in line with 
the log divergent Ximp^ f° r t nc single vacancy. 

The very sudden crossover from divergent to almost 
T independent behavior seen in Fig. llOf al speaks for a 
component Xj^p aligning strongly to the local Neel or- 
der (which becomes the global order at the L depen- 
dent crossover temperature), and justifies the separa- 
tion into parallel and transverse (with respect to the lo- 
cal fluctuating Neel vector) impurity susceptibility com- 
ponents already at intermediate T. However, the lon- 
gitudinal component is not strictly Sf/3T; the recent 
field theory by Sachdev and Vojta predicts that the 
remaining longitudinal contributions, once this leading 
term has been subtracted, has a temperature dependence 
cx Tln(l/T). Nevertheless, the transverse contribution, 
which is cx In (1/T) at low T, dominates. 

In Fig. llUf b). SSE data is shown for the case of the 
two vacancies being nearest neighbors, r = (1, 0), as well 
as at separation r = (2, 1) on the square lattice. Again, 
the single-vacancy data are reproduced at high T. In 
contrast to the divergent trend seen in the maximum- 
separation data in Fig. llUf b) , the finite-size behavior has 
now converged to a near constant at low T, and no signs 
of a log divergence as a function of L is observed. In the 
figure we show only L — 32 results, which are almost con- 
verged to the thermodynamic limit. The absence of log 
corrections for two vacancies at fixed separation is con- 
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FIG. 11: Impurity susceptibilities for different system sizes 
L of the models with a missing bond (solid symbols) and 
a ferromagnetic bond with Jf — J (open symbols). The 
dotted curve shows the expected high-T behavior 1/6T for 
large Jf/J- 

sistent with results of a Green's function calculation, 16 
where the introduction of a second extra spin destroyed 
the log divergence in the frequency dependent T = sus- 
ceptibility observed for the system with a single extra 
spin. 



D. Si = impurities 

We next turn to the QMC results shown in Fig. II II for 
the 2D Heisenberg antiferromagnet with a ferromagnetic 
bond or a missing bond, i.e., with Jp = J or Jp = 0, 
respectively. In this case the impurity susceptibilities do 
not diverge as T — > 0, and the results are not, there- 
fore, multiplied with T. The observed high-T behavior 
of each model, xfmp ~ * 0' 1S due to the fact that the 
susceptibilities of the pure and the doped models can- 
cel, since there is an equal number of independent spins 
in both models at high temperatures. For the missing- 
bond impurity, low-T finite-size effects are clearly seen 
for L = 4 and 8, while for the largest system size L = 32, 
the results should be almost size-converged and show 
little temperature dependence at low T. The observed 
finite-size behavior, xf mp (T — > 0) — > reflects the S = 
ground state, and clearly the size-converged T depen- 
dence also speaks for an Si = impurity. Results for 
the ferromagnetic-bond impurity arc limited to temper- 
atures down to T/J = 1/8, because of the sign problem 
caused by the frustrating ferromagnetic bond. The data 
are reminiscent of the missing-bond results, and hence 
also the ferromagnetic-bond impurity has Si = 0. Both 
models are, clearly, special cases of the system with one 



ferromagnetic bond of arbitrary strength Jp. It would 
be interesting to investigate how the impurity spin mag- 
nitude Si changes as Jp is increased. For Jp/J ^> 1, 
the two spins connected by the ferromagnetic bond form 
a triplet and hence should give an Si = 1 Curie con- 
tribution Si(Si + 1)/3T = 2/3T when J < T < J F . 
The remaining N — 2 spins each contribute 1/4T, and 
hence the impurity susceptibility should be 1/6T in this 
regime. In Fig.^Jthe results for T > J are closer to this 
form for Jp — J than for Jp — 0, but the requirement 
J < T < Jp is not satisfied and the deviations (reduction 
relative to 1/6T) reflect an expected crossover from the 
high-T independent-spin form xf mp ~ 0. 

An interesting question is whether the classical-like 
Curie behavior xf mp ~ Sf/3T with Si = 1 can be ob- 
served in this model for Jp > J. As already discussed in 
Sec. II, the asymmetric coupling to the bulk of the two 
spins connected by the ferro bond most likely implies a 
T — * behavior corresponding to Si = for any finite 
Jp. This is because an Si — 1 impurity requires that 
the two spins at the ferro bond are dominantly in the 
m z = 1 state | |T) with respect to the local Neel order 
(in a semiclassical picture such as our effective impurity 
model), whereas in fact the couplings in this case instead 
favor the m z = component (| U) + I IT>)/\/2— 

E. Internal energy and specific heat 

We finally discuss our SSE calculations concerning im- 
purity effects on the internal energy and the specific heat, 
which we have obtained using the estimators in Eqs. I|l(ja|) 
and ijl6b(> . respectively. In analogy to Eq. J2J), we again 
define the impurity quantities as differences between the 
doped and the pure systems, i. e., 

El2 P = E { i ) -E (a) , (23a) 

C Sp = C {i) ~ C (a)i ( 23b ) 

where i — b,c, d, e, and /, correspond to the different 
impurity systems shown in Figs.nand[21 and the symbol 
a denotes the pure system. In Fig. I12l results for the im- 
purity energies are shown for the vacancy model (a) and 
the added-spin model (b) with J± = J and J± = J/2 
(shown in the inset). At high T the impurity energies 
vanish, since the mean energy of each independent spin 
becomes zero. For L — 64, all finite-size effects are elim- 
inated within statistical errors for both models in the 
T range considered. Since the vacancy system has four 
antiferromagnetic bonds less than the pure system, the 
impurity energy E^ p , shown in Fig. I12f a). is positive 
at all T. At low T the results converge to a constant 
value, which should be equal to the energy cost of re- 
moving one spin from an infinite lattice in its ground 
state. The low-T value observed in Fig. I12f a) is indeed 
consistent with T = results obtained in a previous lin- 
ear spin-wave study^ Results for the added-spin model 
with J± = J, shown in Fig. I12f b). are negative because 
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FIG. 12: The impurity energies of the (a) vacancy and (b) 
added-spin (Jj_ = J) models, for different system sizes L. The 
inset shows the L = 64 QMC results for the added-spin model 
with J± — J/2. The curves are results of the corresponding 
effective models. 



FIG. 13: The impurity specific heats of the (a) vacancy and 
(b) added-spin (Jj_ = J) models, for different system sizes 
L. The dashed curve in (a) shows the result of the effective 
model. 



of the one extra antiferromagnetic bond, and the size- 
converged behavior seems to also tend to a constant as 
T — > 0. This constant value corresponds to the energy 
cost of removing the off-plane added spin from its host 
lattice, and its magnitude is observed to be roughly one 
fourth of the low-T value of the vacancy impurity en- 
ergy. The T dependence of the L — 64 results for the 
added-spin impurity with J± = J/2, shown in the inset, 
are qualitatively very similar, but because of the smaller 
impurity-bond strength the absolute values arc smaller. 

The solid curves in Fig. ^| are results of the corre- 
sponding effective models. Using the same values on a/ J 
and r/J as previously when calculating the impurity sus- 
ceptibilities, we obtain a qualitative agreement for the 
vacancy model while the agreement is remarkably good 
for the added-spin model, both for J± = J and J± = J/2 
(inset). Hence, in addition to reproducing the impurity 
susceptibilities of the full models, the effective models 
also describe properly the energetics of the full models. 
Also, the parameters a and r can be tuned to give a bet- 
ter agreement for the vacancy model in Fig. I12f al. but 
this in turn will give a poorer agreement between QMC 
and effective-model results for the impurity susceptibility 
of the vacancy model in Fig. |SJ 



In Fig. El QMC results for the impurity specific heats 
are shown for the vacancy model (a) and the added- 
spin model (b) with J± = J. As the system size is 
increased and the temperature is lowered the statistical 
errors grow rapidly. The size-converged behavior is diffi- 
cult to determine below T/J w 0.3, but C^ p in Fig.GJa) 
is, nevertheless, consistent with the behavior of in 

Fig. da), as C = dE(T)/dT. The point at which C^ p 
goes trough zero, T/J « 0.5, corresponds to the maxi- 
mum in the energy curve in Fig. ^| The effective 
model reproduces well the high-T behavior and also ex- 
hibits a negative minimum at intermediate temperature. 
However, this features is much less pronounced than for 
the full model, and the maximum at lower T is missing. 
In Fig. I12f b). sufficient accuracy in the simulations has 
not been reached for larger system sizes, and the size- 
converged behavior can therefore not be determined. 



V. SUMMARY 

In this paper we have presented results of an extensive 
QMC study of impurity effects in the 5=1/2 Heisen- 
berg antiferromagnet on a square lattice, as well as some 



14 



results for a 3D system. The effects of different types 
of single static impurities on the magnetic susceptibility 
and the specific heat have been investigated. 

For several types of Si ^ impurities in 2D (vacancy, 
added spin, ferromagnetically coupled spin), our very 
precise simulation data has revealed an additive loga- 
rithmic correction to the predicted classical-like Curie 
contribution Sf/3T to the impurity susceptibility. We 
have argued that this logarithmic contribution reflects 
primarily fluctuations transverse to the local Neel order 
at the impurity. This is in agreement with recent field- 
theoretical workjS&Sl carried out after our initial report 
of log corrections^ Here we have shown that our numer- 
ical results are in excellent quantitative agreement with 
these field-theoretical results?^^ containing both lead- 
ing and subleading logarithmic corrections. In 3D, we 
find no signs of logarithmic corrections, in accord with 
predictions^ 

In order to have a simple mechanism explaining the 
leading-order (i.e., apart from the log corrections) impu- 
rity physics, we have also introduced few-spin effective 
models. Comparisons with the QMC results show that 
the effective models can distinguish between impurities 
of different types and spins Si. In many cases the quan- 
titative agreement between the effective and full models 
is surprisingly good over a wide temperature range. This 
suggests that extended effective models based on larger 
clusters of spins, e.g., 3x3 clusters centered around the 
site impurities, should give very accurate descriptions, 
perhaps also for the vacancy model which we here found 
was the hardest case to describe with the simplest effec- 
tive model. 
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APPENDIX: SELECTED QMC DATA FOR THE 
SUSCEPTIBILITY, ENERGY, AND SPECIFIC 
HEAT 

The numerical data underlying the analyses carried out 
in this paper are of very high accuracy — the small errors 
are only statistical in nature — and may hence be useful as 
bench marks for alternative calculations. In Tables [I] and 



[n]we therefore list L — 64 data for the susceptibility and 
the internal energy, at several inverse temperatures J/T, 
for the pure (a), vacancy (b), and added-spin models (c). 

In Fig. E| we show the SSE results for the specific heat 
of the pure 2D Heisenberg antiferromagnet at tempera- 
tures down to Tj J = 1/32. At such low temperatures the 
specific heat has not been determined reliably in previous 
studies. 43 We have obtained the results using the direct 
estimator, Eq. I|16b(l . The low-T data shown in the in- 

TABLE I: Selected L = 64 data for X(i)/L 2 at inverse tem- 
perature J/T, where i — a,b, and c correspond to the pure, 
vacancy, and added-spin systems, respectively. 

i — c 

J/T i = a i = b J± = J J± = J/2 

32 0.045043(3) 0.045767(4) 0.045907(4) 0.046128(3) 

16 0.046359(3) 0.046737(3) 0.046867(3) 0.047056(3) 

8 0.049144(3) 0.049344(3) 0.049461(2) 0.049557(3) 

4 0.055994(1) 0.056092(1) 0.056179(1) 0.056214(1) 

2 0.0786001(4) 0.0786347(4) 0.0786945(4) 0.0787106(4) 

1 0.0935393(2) 0.0935377(2) 0.0935859(2) 0.0935937(2) 



TABLE II: Selected L = 64 data for -E (i) /L 2 at inverse 
temperature J/T, where i — a,b, and c correspond to the 
pure, vacancy, and added-spin systems, respectively. 

i — c 

J/T i = a j = b Jx = J J± = J/2 

32 0.6694416(5) 0.6691562(5) 0.6695154(6) 0.6694702(5) 

16 0.6693890(7) 0.6691048(7) 0.6694625(7) 0.6694158(7) 

8 0.6689102(9) 0.6686247(8) 0.668983(1) 0.6689330(9) 

4 0.663745(1) 0.663452(1) 0.663807(1) 0.663761(1) 

2 0.593051(1) 0.592755(1) 0.593093(1) 0.593060(2) 

1 0.387560(2) 0.387372(2) 0.387600(2) 0.387569(2) 



set of Fig. E] are clearly consistent with the quadratic T 
behavior suggested in the Hasenfratz-Niedermeyer chiral 
perturbation theory i£i 

C(T) = ^$-T 2 + 0(T 4 ), (A.l) 

where we use c = 1.66 for the spin- wave velocity 4 ^ and 
C(3) = 1.2020569. 
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FIG. 14: Size-normalized specific heats of the pure 2D Heisen- 
berg antiferromagnet for different system sizes L. Error bars 
are smaller than the symbols. The inset compares our low-T 
data with theory (solid curve) ( Ref 4i). 
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